An increased estimate of the merger rate of double neutron 
stars from observations of a highly relativistic system 
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7— I ■ 

■ The merger 1 of close binary systems containing two neutron stars should produce a burst of gravi- 

tational waves, as predicted by the theory of general relativity 2 . A reliable estimate of the double- 
neutron-star merger rate in the Galaxy is crucial in order to predict whether current gravity wave 
detectors will be successful in detecting such bursts. Present estimates of this rate are rather low 3-7 , 
because we know of only a few double-neutron-star binaries with merger times less than the age of 
the Universe. Here we report the discovery of a 22-ms pulsar, PSR J0737 3039, which is a member 
of a highly relativistic double-neutron-star binary with an orbital period of 2.4 hours. This system 
will merge in about 85 Myr, a time much shorter than for any other known neutron-star binary. To- 
gether with the relatively low radio luminosity of PSR J0737 3039, this timescale implies an order-of- 
magnitude increase in the predicted merger rate for double-neutron-star systems in our Galaxy (and 
» ■ in the rest of the Universe) . 

c§ : 

PSR J0737— 3039 was discovered during a pul- band gave an improved position and flux density 
sar search carried out using a multibeam receiver 8 for the pulsar. Knowledge of the pulsar position 
on the Parkes 64-m radio telescope in New South with subarcsecond precision allowed determination 
Whales, Australia. The original detection showed a of the rotational period derivative, P, and other pa- 
large change in apparent pulsar period during the 4- rameters from the available data span. Table 1 re- 
min observation time, suggesting that the pulsar is a ports results derived from a coherent phase fit to 
member of a tight binary system. Follow-up obser- data taken over about five months. The measured 
vations undertaken at Parkes consisting of continu- value of Co = 16.88 yr _1 is about four times that 
ous ~ 5-hour observations showed that the orbit has of PSR B1913+16 (ref. 9), previously the highest- 
a very short period (2.4 hrs) and a significant ec- known. If the observed Co is entirely due to gen- 
centricity (0.088). The derived orbital parameters eral relativity, it indicates a total system mass M — 
implied that the system is relatively massive, proba- 2.58 ± 0.02 M , where M Q is the mass of the Sun. 
bly consisting of two neutron stars, and predicted a Figure 1 shows the constraints on the masses of 
huge rate of periastron advance Co due to effects of the pulsar and its companion resulting from the ob- 
general relativity. Indeed, after only a few days of servations so far and the mean pulse profile as an 
pulse-timing observations we were able to detect a inset. The shaded region indicates values that are 
significant value of Co. ruled out by the mass function Mf and the observed 

Interferometric observations made using the Aus- Co constrains the system to lie between the two di- 

tralia Telescope Compact Array (ATCA) in the 20-cm agonal lines. Together, these constraints imply that 
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the pulsar mass m p is less than 1.35 M Q and that 
the companion mass m c is greater than 1.24 M Q . 
The derived upper limit on m p is consistent with the 
mass distribution of pulsars 10 in double-neutron-star 
systems. The minimum allowed companion mass 
implies that the companion is very likely also to be a 
neutron star. 

This is supported by the observed significant ec- 
centricity of the orbit, indicative of sudden and/or 
asymmetric mass loss from the system in the ex- 
plosion when the companion collapsed 11 . Sig- 
nificant classical contributions to lu, either from 
tidal deformation 12 or rotation-induced quadrupole 
moment 13 ' 14 in the companion - effects which might 
be present if it were not a neutron star - are there- 
fore unlikely. In fact, lj gr <c lo would imply an un- 
acceptably low value for the pulsar mass. Hence this 
system probably contains two neutron stars with rel- 
atively low masses in the range 1.24 - 1.35 M Q , and 
the binary orbit is probably nearly edge-on (inclina- 
tion i w 90°). For i = 68°, the implied mass of the 
pulsar is 1.24 M Q and that of the companion is 1.35 
M Q . 

A precise determination of the masses of both 
stars, and then of the orbital inclination, will be ob- 
tained when the gravitational redshift 9 , 7, is mea- 
sured. The expected minimum value (assuming an 
edge-on orbit) of this parameter is 383 fj,s. How- 
ever, its orbital dependence is covariant with a sin i 
(Table 1) until there has been significant periastron 
precession. Whereas we currently get a formal value 
of 7 = 400 ± 100 [is from a fit of the timing data, a 
useful measurement of 7 should be possible within a 
year. 

The parameters of the PSR J0737-3039 binary 
system, and the high timing precision made possible 
by the narrow pulses of this pulsar, promise to make 
this system an excellent 'laboratory' for the investi- 
gation of relativistic astrophysics, in particular of all 
the effects that are enhanced by a short orbital pe- 
riod Pb- Geodetic spin-precession due to curvature 
of space-time around the pulsar is expected to have 
a period of 75 yr. This is four times less than the pre- 
viously known shortest spin precession period and 
should allow detailed mapping of the pulsar emis- 
sion beam 15 in few years. It may also have the unde- 
sirable effect of making the pulsar undetectable in a 
short time. 

The secular decay of the orbital period due to 
gravitational wave emission should be measurable 



within about 15 months. In PSR B1913+16, uncer- 
tainty in the orbital period derivative due to acceler- 
ation in the Galactic potential limits the precision of 
this test of the predictions of general relativity 1617 . 
By comparison, PSR J0737— 3039 is relatively close 
to the Sun and the uncertainty in Galactic accelera- 
tion will be much smaller, allowing a more precise 
test. Because of the probable high orbital inclina- 
tion of the system, the delay due to the curvature of 
the space-time in the vicinity of the companion (the 
Shapiro delay 18 ) should be measurable. 

The relatively small orbital eccentricity suggests 
that this system may have already undergone a sub- 
stantial change of the orbital elements due to grav- 
itational decay. Assuming that the evolution of or- 
bital separation and eccentricity are entirely driven 
by gravitational radiation emission 19 , and adopting 
a value of 100 Myr for the time elapsed since birth 4 , 
Tb (see Table 1), the orbital initially had a period of 
3.3 h and an eccentricity of 0.119. Hence this sys- 
tem would already have experienced a ~ 25% de- 
cay of the orbital period and circularization, which 
would be the largest variation of these parameters 
in the sample of the known double-neutron-star bi- 
nary systems. 

Of the five double-neutron-star systems known 
so far 20 , only three have orbits tight enough that the 
two neutron stars will merge within a Hubble time. 
Two of them (PSR B1913+16 and PSR B1534+12) 
are located in the Galactic field, while the third 
(PSR B2 12 7+ 11C) is found on the outskirts of a 
globular cluster. The contribution of globular cluster 
systems to the Galactic merger rate is estimated to 
be negligible 21 . Also, recent studies 6 have demon- 
strated that the current estimate of the Galactic 
merger rate 1Z relies mostly on PSR B1913+16. So, 
in the following discussion, we start by compar- 
ing the observed properties of PSR B1913+16 and 
PSR J0737-3039. 

PSR J0737-3039 and its companion star will 
merge due to the emission of gravitational waves in 
T m s=a 85 Myr, a timescale that is a factor 3.5 shorter 
than that for PSR B1913+16 (ref. 9). In addition, 
the estimated distance for PSR J0737-3039 (500 - 
600 pc, based on the observed dispersion measure 
and a model for the distribution of ionised gas in the 
interstellar medium 22 ) is an order of magnitude less 
than that of PSR B1913+16. These properties have 
a substantial effect on the prediction of the rate of 
merging events in the Galaxy. 
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For a given class k of binary pulsars in the 
Galaxy, apart from a beaming correction factor, the 
merger rate TZ^ is calculated 6 as 1Z k oc N k /r k . 
Here r k is the binary pulsar lifetime defined as 
the sum of the time since birth 4 , Tb, and the re- 
maining time before coalescence, r m , and N k is 
the scaling factor defined as the number of bina- 
ries in the Galaxy belonging to the given class. The 
shorter lifetime of PSR J0737-3039 (V1913/T0737 = 
(365 Myr)/(185 Myr) « 2, where the subscript num- 
bers refer to the pulsars), implies a doubling of 
the ratio 7?.o737/7?.i9i3- A much more substantial 
increase results from the computation of the ratio 
of the scaling factors ^0737/^1913. The luminosity 
L 400 w 30 mJy kpc 2 of PSR J0737-3039 is much 
lower than that of PSRB1913+16 (~ 200mJykpc 2 ). 
For a planar homogeneous distribution of pulsars in 
the Galaxy, the ratio N 0737 /N 1913 oc £1913/^0737 ~ 
6. Hence we obtain ^0737/^1913 ~ 12. Includ- 
ing the moderate contribution of the longer-lived 
PSR B1534+12 system to the total rate 4 " 7 , we ob- 
tain an increase factor for the total merger rate 

C&0737 + ^1913 + ftl534)/(ftl913 + ^1534) Of about an 

order of magnitude. A better estimate of this factor 
and its uncertainty can be obtained using a bayesian 
statistical approach and accounting for the full lumi- 
nosity function of pulsars in the Galaxy 6 . 

Figure 2 dispalys an estimate of the probability 
density function for the merger rate increase factor, 
showing a peak value of ~ 8 and an upper limit of 
~ 30 at a 95% confidence level. This result, derived 
in the simple assumption of a fixed pulsar luminos- 
ity and a uniform disk distribution of pulsar bina- 
ries, can be refined by including the parameters of 
the present survey into a simulation program mod- 
elling survey selection effects (and hence detection 
probability) and the Galactic population of pulsars. 
The Parkes high-latitude survey covers the region 
enclosed by Galactic latitude \b\ < 60° and Galac- 
tic longitude 220° < I < 260° for a total of 6456 
pointings, each of duration 265 s, using a band- 
width of 288 MHz and a sampling time of 125 (is. 
PSR J0737-3039 was detected with signal-to-noise 
ratio of 18.7 in a standard Fourier search 23 . Varia- 
tion in the apparent pulsar period due to binary mo- 
tion during the discovery observation resulted in a 
30% reduction in the observed signal-to-noise ratio. 

Extensive simulations (VK. et ah, manuscrtpt in 
preparation) produce results consistent with that de- 
rived in Fig. fig2, and show that the peak of the 



merger rate increase factor resulting from the dis- 
covery of PSR J0737— 3039 lies in the range 6 to 
7 and is largely independent of the adopted pul- 
sar population model. On the other hand, the ac- 
tual predicted value of the Galactic merger rate and 
hence the detection rate by gravity wave detectors 
depends on the shape of the pulsar luminosity func- 
tion. For the most favourable distribution model 
available (model 15 of ref. 6), the updated cosmic 
detection rate for first-generation gravity wave de- 
tectors such as VIRGO 24 , LIGO 25 and GEO 26 can be 
as high as 1 every 1-2 years at 95% confidence level. 

After a few years of operation of the gravity wave 
detectors, it should be possible to test these predic- 
tions directly and thus place better constraints on the 
cosmic population of double-neutron-star binaries. 

Received 12 August; accepted 15 October 2003. 
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1.5 

Pulsar Mass (M ) 

Figure 1 Constraints on the masses of the neutron star and the 
companion. The keplerian orbital parameters give the mass func- 
tion: 



M f 



(m c sin i) 
(m c + ra p ) 2 



4-7T 2 (a sini) 3 
GP 2 b 



O.291O54(8)M 



where m p and m c are the masses of the pulsar and companion, 
Pt, is the orbital period, and a sini is the projected semi-major 
axis of the pulsar orbit; the estimated uncertainty in the last 
quoted digit is given in parentheses. This provides a constraint on 
the minimum m c as a function of rrt p by setting i = 90°, where 
i is the inclination of the orbit normal to the line of sight. In 
addition, assuming that the periastron advance rate ui is entirely 
due to general relativity, (that is, u> = ugr), its value can be 
combined with the other orbital parameters to determine the total 
system mass M = m c + m p using the equation 27 : 



p 5/2 
f_b 

2-kG 



(l-e 2 )c 2 ^ 



-i 3/2 



6tt 



2.58(2)M 



where e is the orbital eccentricity. In the m c — m p plane, the 
shaded region indicates the mass values ruled out by the mass 
function Mf, and the intersection of the boundary of this region 
with the ill limits (defining a diagonal strip in the plot) provides a 
further constraint on the masses. The dash-dotted lines represent 
the mass values corresponding to the relativistic parameter 7 for 
an edge-on orbit (7 a ) and for an orbital inclination i = 68° (7b)- 
The inset shows the mean pulse profile at 1400 MHz from a 
5-h observation. There are two pulse components separated by 
approximately half the 22-ms period. 




1 10 
Merger rate increase factor 

Figure 2 Probability density function for the increase in 
the double-neutron-star merger rate (72.0737 + 7?.igi3 + 
7£i534)/(7£i9i3 + 721534) resulting from the discovery of 
PSR J0737— 3039. For a given class k of binary pulsars in the 
Galaxy, the probability density function P(72k) for the corre- 
sponding merger rate TZ^ is obtained from the relation 6 P(72. k ) = 
A 2 Tl} i e~ A ' R - k where A = Tk/(AT k / k ), AT k is a population scaling 
factor (see text), r k is the binary pulsar lifetime and / k is the cor- 
rection factor due to the beamed nature of the radio pulsar emis- 
sion. The dashed line represents the reference probability density 
function corresponding to the merger rate 72 id = 72i9i3 +7^1534 
due to PSR B1913+16 and PSR B1534+12 only 6 . The heavy 
solid line represents the probability density function correspond- 
ing to the new merger rate 72-ncw = 72o737 + 72-1913 + 72.1534. 
The parameters adopted for PSR J0737— 3039 are those derived 
in the text (-/V0737 = 6./V1913, T0737 = 185 Myr) and the beam- 
ing factor for PSR J0737— 3039, /0737, is chosen as the average 
of the beaming factor of the other two known merging binaries 
(/0737 = l-06/i9i3). The dotted vertical lines represent the 68% 
and 95% confidence level boundaries on the determination of the 
increase factor. The dominant role of PSR J0737— 3039 in shap- 
ing the new statistics of the double-neutron-star merger rate is 
evident. 
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6.9(6) 


Flux density at 430 MHz (mjy) 


100(20) 


Period, P (ms) 


22.69937854062(6) 


Period derivative, P (s/s) 


2.3(6) x 10 18 


Epoch (MJD) 


52800.0 


Orbital period, P b (days) 


0.102251561(8) 


Projected semi major axis, asini (It-s) 


1.415176(5) 


Eccentricity 


0.087793(8) 


Longitude of periastron u> 


68°. 743(9) 


Epoch of periastron, TO (MJD) 


52760.807391(3) 


Advance of periastron ill (°yr _1 ) 


16.88(10) 


Post-fit residual (/is) 


33 


Mass function, M f (M Q ) 


0.291054(8) 


Total system mass, M (Mq) 


2.58(2) 


Characteristic age, r c (Myr) 


160(50) 


Age since birth Tb (Myr) 


100(50) 


Surface dipole magnetic field strength, B (G) 


7.3(9) x 10 9 


Expected gravitational redshift parameter, 7 (fis) 


> 383 


Expected orbital period derivative (s s _1 ) 


-1.24 x 10" 12 



Table 1 Measured and derived parameters for PSR J0737 3039. The 
pulsar position is obtained from observations made with the ATCA used 
in pulsar-gating mode. Uncertainties are given in parentheses and refer 
to the last quoted digit. The flux at 1400 MHz is derived from the ATCA 
observations, while the flux at 430 MHz is obtained from Parkes obser- 
vations with the flux value being derived from the observed signal-to- 
noise ratio and the estimated receiver sensitivity. Best-fit timing param- 
eters are derived from arrival-time data by minimising the x 2 function 
using the program TEMPO (see http://pulsar.princeton.edu/tempo) 
with the relativistic binary model of ref. 27. For all the timing parame- 
ters except Pb, u> and P, uncertainties are twice the formal-fit la error. 
Because of the short time span, u> is strongly covariant with other orbital 
parameters and the nominal la error resulting from the phase coherent 
fit of the pulse time-of-arrival might not represent the actual lcr sta- 
tistical uncertainty. By producing contours of constant \ 2 values as a 
function of Pj, and u>, we have verified that an uncertainty of the order 
of four times the nominal lcr error represents the uncertainty range for 
these parameters. The uncertainty in P is due to the residual covari- 
ance with the position terms resulting from the available uncertainty 
in the celestial coordinates obtained at the ATCA. The characteristic 
age r c and the surface dipole magnetic field strength B are derived 
from the observed P and P using the relations 28 r c = 0.5P/P and 
B = 3.2 X 10 19 (PP) 1 / 2 . The age since birth, r b , is the time since the 
binary pulsar left the spin-up line 4 . 



